
from numpy import *
from pylab import *
import Well
import ExperimentalAnalyzer as analyze
import PIL
from PIL import Image
import os

""" This class represents well maps, providing a means to build and manipulate them, as well as get
relevant data from specific instances. """
class wellMap:
    WELL = 1

    def __init__(self, plateDimensions, bMatrix, gThreshold=2.0, rThreshold=-1.0):
        # Notice that rThreshold is defaulted only so you can run the program without specifying. If you want do
        # do anything useful involving talkers, you need to adjust this.
        self.wellMatrix=zeros(plateDimensions, 'int')
        self.wellList=[] # Holds set of well coordinates.
        self.tNot=[] # This holds the intensities - in a one-to-one mapping with wellList -- of each
                     # each well at the initial image. Used in calculating the average background intensity
                     # of images.
        self.numWells = 0
        self.bgMatrix = bMatrix
        self.greenThreshold=gThreshold 
        self.redThreshold=rThreshold
        self.redWells=[]
        self.greenWells=[]
        self.sytosImg = None # Only one 'red' - sytos AKA talkers - image per experiment. Set this before doing anything with reds.
        print "Well Map initiated with plateDimensions: (" + str(plateDimensions[0]) + ", " + str(plateDimensions[1]) + ")"
        
    def setSytosImg(self, img):
        """ Setter for this experiment's associated 'talker' image. """
        self.sytosImg = img
        
    def findAllTalkerWellsWithF(self, f, concentration, slack=200):
        print "Finding talker wells using estimation from concentration.."
        missedVal, initially = False, ""
        curThresh = 1500 # Totally arbitrary starting threshold
        adjustFactor = 300 # Also somewhat arbitrary
        metaAdjustFactor = (1.3, .7) # Adjustment scalar for the adjust factor
        if len(self.wellList)== 0:
            print "No wells to consider. Build the wellList before calling."
            return
        if self.sytosImg == None:
            print "Sytos image is not set! Please set and try to run again."
            return
        expectedRedCount = f(concentration)
        print "Expected red count is %s." % str(expectedRedCount)
        while (abs(len(self.redWells) - expectedRedCount) > slack):
            print "Reds found so far: %d" % len(self.redWells)
            print "Of by: %s" % str(abs(len(self.redWells) - expectedRedCount))
            print "Slack: %s" % str(slack)
            print "Adjust factor is %s" % str(adjustFactor)
            #print "MetaAdjust factor is %s" % str(metaAdjustFactor)
            print "curThresh is currently %s" % str(curThresh)
            # While the absolute difference between the expected count and
            # the found count is greater than slack, adjust the threshold
            self.findAllTalkerWells(rThresh=curThresh)
            if initially=="":
                if self.redWells < expectedRedCount: 
                    initially = "less than expected"
                else:
                    initially = "more than expected"
            print "Initially is: " + initially
            # Now adjust the adjusting var
            if len(self.redWells) < expectedRedCount:            
                # If we were too low, then the threshold was too high
                # if we initially found more than we expected and
                # now we're finding less, we've missed the value
                if initially == "more than expected":
                    # If we miss the value, lower our adjustFactor and
                    # set the initially var to less than; so that
                    # we can re-adjust if we go past it again
                    print "We missed it! Lowering adjustFactor"
                    adjustFactor = .75*adjustFactor
                    initially = "less than expected"
                curThresh = curThresh - adjustFactor
                print "Too Low! curThresh adjusted to: %s" % str(curThresh)
            elif len(self.redWells) > expectedRedCount:
                # Then we were too high and need to increase the threshold
                # If we initially found less than expected and now we're finding
                # more, we've missed the value.
                if initially == "less than expected":
                    print "We missed it! Lowering adjustFactor.."
                    adjustFactor = .75*adjustFactor
                    initially = "more than expected"
                curThresh = curThresh + adjustFactor
                print "Too High! curThresh adjusted to: %s" % str(curThresh)
            
            """
            # Here we adjust the adjust factor
            if abs(len(self.redWells) - expectedRedCount) > 2*adjustFactor:
                # We're way off, increase the adjust factor
                adjustFactor = adjustFactor * metaAdjustFactor[0]
                print "Way Off!! Meta-Adjusting.. adjustFactor is now %s" % str(adjustFactor)
            if abs(len(self.redWells) - expectedRedCount) < adjustFactor:
                # We're pretty close, decrease the adjust factor
                adjustFactor = adjustFactor * metaAdjustFactor[1]
                print "Pretty Close... Meta-Adjusting.. adjustFactor is now %s" % str(adjustFactor)
            """
            print "\n"
        
    def findAllTalkerWells(self, rThresh=-1, boundsBox=None, discountConflictWells=True):
        """ Populates the redWellLs list with red wells. 
        
        Requisite parameters:
        n -- If a well has this many 'red' pixels (a pixel is read if it exceeds self.redThreshold
        
        Optional parameters:
        discountConflictWells -- If this is True (as is by default) then a well will be ignored if it contains both a red and a green pixel.""" 
        print "findAllTalkerWells -- rThresh is %s"%str(rThresh)
        print "Clearing red well list.. "
        self.redWells =[]
        
        if rThresh == -1:
            rThresh = self.redThreshold    
            
        if len(self.wellList)== 0:
            print "No wells to consider. Build the wellList before calling."
            return
        
        if self.sytosImg == None:
            print "Sytos image is not set! Please set and try to run again."
            return
        
        for well in self.wellList:
            redIntense = self.getAvgIntensityForWell("red", well, self.sytosImg)
            #gIntens = self.getAvgIntensityForWell("green", well, img)
            if redIntense > rThresh:
                self.redWells.append(well)
        print "Found " + str(len(self.redWells)) + " red wells"
        if not boundsBox == None:
            print "Filtering red wells: returning only those within " + str(boundsBox)
            self.redWells = filter(lambda well: analyze.pointFallsWithinBounds((well.x,well.y), boundsBox), self.redWells)
        return self.redWells
        
    
    # Need to have a cut off bounds
    def findAllGreensUsingCutOff(self, imgDir, imgCutOff, imgStart=0, alreadyResponderLs=[], boundsBox=None,cutOffFactorBounds=(0,900)):
        print "Finding all green wells in well map (using CUT-OFF method)..."
        candidates = filter(lambda x: x not in self.redWells, self.wellList)
        if not boundsBox == None:
            print "Bounding search region for green wells."
            candidates = filter(lambda w: analyze.pointFallsWithinBounds((w.x, w.y), boundsBox), candidates)
        print "Considering " + str(len(candidates)) + " candidate green wells..."
        self.greenWells = []
        imgList=os.listdir(imgDir)
        imgList.sort()
        for image in imgList[imgStart:imgCutOff]:
            i = None
            print "Opened: " + image
            i = Image.open(imgDir + "/" + image)
            for w in candidates:
                avgIntens = self.getAvgIntensityForWell("green", w, i)
                #print avgIntens
                """
                print "avg intensity: " + str(avgIntens)
                print cutOffFactorBounds[0]
                print cutOffFactorBounds[1]
                """
                if avgIntens >= cutOffFactorBounds[0] and avgIntens <= cutOffFactorBounds[1] and not w in alreadyResponderLs:
                    #if not w in self.greenWells:
                    self.greenWells.append(w)
                    print "A well in the following image has exceeded or met the cutoff: " + image
            #self.greenWells = set(self.greenWells)
        greenWellsP = []
        for gW in set(self.greenWells):
            greenWellsP.append(gW)
        self.greenWells = greenWellsP
        #return set(self.greenWells)
        return self.greenWells
        
        
    def findAllGreensUsingStdDev(self, imgDir, m, imgCutOff=50, boundsBox=None):
        print "Finding all green wells in map (using standard deviation method)"
        candidates = filter(lambda x: x not in self.redWells, self.wellList)
        if not boundsBox == None:
            print "Bounding search region for green wells."
            candidates = filter(lambda w: analyze.pointFallsWithinBounds((w.x, w.y), boundsBox), candidates)
        print "Considering " + str(len(candidates)) + " candidate green wells..."
        self.greenWells = []
        for w in candidates:
            continue
            
        
    def findAllGreens(self, imgDir, imgCutOff=50, boundsBox=None):
        # Takes base-delta to be the max of the k-lowest deltas.
        k = 50
        print "Finding all green wells in well map (using DELTA method)..."
        d={}
        candidates = filter(lambda x: x not in self.redWells, self.wellList)
        if not boundsBox == None:
            print "Bounding search region for green wells."
            candidates = filter(lambda w: analyze.pointFallsWithinBounds((w.x, w.y), boundsBox), candidates)
        print "Considering " + str(len(candidates)) + " candidate green wells..."
        for w in candidates:
            d[str(w.x)+","+str(w.y)]=[]
        print "Constructed dictionary."
        imgList=os.listdir(imgDir)
        for image in imgList[0:imgCutOff]:
            i = None
            print "Opened: " + image
            i = Image.open(imgDir + "/" + image)
            for w in candidates:
                avgIntens = self.getAvgIntensityForWell("green", w, i)
                d[str(w.x)+","+str(w.y)].append(avgIntens)
   
        kLowestDelta = []
        for wKey in d.keys():
            intensityLs = d[wKey]
            delta = max(intensityLs) - min(intensityLs)
            #print delta
            d[wKey]=delta
            if len(kLowestDelta)<k:
                #kLowestDelta.append([wKey, delta])
                kLowestDelta.append(delta)
            elif delta < min(kLowestDelta):
                kLowestDelta[kLowestDelta.index(min(kLowestDelta))] = delta
            """
            else:
                kDeltaLs = map(lambda x:x[1], kLowestDelta)
                if delta < min(kDeltaLs):
                   kLowestDelta[kDeltaLs.index(min(kDeltaLs))]=([wKey, delta])
           """
        print kLowestDelta
        deltaBase = max(kLowestDelta)
        #deltaBase = max(map(lambda x:x[1], kLowestDelta))
        # Now iterate through the list; we'll consider a well 'green' if it's
        # delta exceeds delta base times some scalar.
        print "DELTA-BASE: " +  str(deltaBase)
        self.greenWells = []
        for w in candidates:
            if d[str(w.x)+","+str(w.y)] > deltaBase:
                self.greenWells.append(w)
        print "Done. Returning a list of " + str(len(self.greenWells)) + " wells considered greens.\n"
        return self.greenWells
    
    def addWell(self, well):
        """ Inserts a well location into the pixel matrix. Note
        that we use the pixels of the plateDimensions as our coordinate
        system. 
        
        Requisite parameters:
        well --- this is the well object that will be added to self.
        """
        self.wellMatrix[well.x][well.y]=self.WELL
        self.wellList.append(well)
        self.numWells += 1
    
    def setTNot(self, img):
        """ Set's this instances tNot image. """
        self.tNotImg = img
         
    def getTNotIntensityForWell(self, w):
        """ Returns the intensity logged for the given well at time t=0 """
        total=0
        for xPrime in (0,1,2):
            for yPrime in (0,1,2):
                try:
                   p = self.tNotImg.getpixel((w.x+xPrime, w.y+yPrime)) 
                except: 
                   # We're out of bounds here, just ignore
                   p = 0
                total += p
        return total/9  
    
    def getKNearestNeighbors(self, well, k):
        """ Returns the locations (upper-left) of the k nearest wells to
        well.
        
        Requisite parameters:
        well --- this is the well object that will be added to self.
        k --- this many neighbor coordinates will be returned.
        """
        nList = []
        count = 0
        # Start sweeping out 1 pixel from origin
        scalar = 1
        y = well.y
        x = well.x
        while count < k:
            count += self.checkPixelRangeK((x,y), scalar, nList)
            scalar += 1
        # In case we got more than k neighbors we just lop off the last ones,
        # which are necassarily further away. 
        return nList[0:k]
    
    def countNumOfGreenWells(self, wellLs):
        """ Returns the number of well objects in wellLs that are also in the green well
        list """
        gCount = 0
        for w in wellLs:
            #print "w " +str(w)
            for rW in self.greenWells:
                #print rW
                if w[0] == rW.x and w[1] == rW.y:
                    gCount += 1
        return gCount
    
    def countNumOfRedWells(self, wellLs):
        """ Returns the number of well objects in wellLs that are also in the red well
        list """
        rCount = 0
        for w in wellLs:
            #print "w " +str(w)
            for rW in self.redWells:
                #print rW
                if w[0] == rW.x and w[1] == rW.y:
                    rCount += 1
        #print "Returning red-count: " + str(rCount)
        return rCount
    
    def checkPixelRangeK(self, p, k, nList):
        """ Check all 8 'surrounding' pixels that are k pixels
        from p for k nearest wells. """
        count = 0
        for xSign in (-1, 0, 1):
            for ySign in (-1, 0, 1):
                # If they were both 0 we'd be counting the origin, p.
                if not xSign==ySign==0:
                    try:
                        xPrime = p[0]+xSign*k
                        yPrime = p[1]+ySign*k
                        if self.wellMatrix[xPrime, yPrime] == self.WELL:
                            nList.append((xPrime, yPrime))
                            count += 1
                    except:
                        continue
        return count
                            
    def getNeighborsInRange(self, xList, yList, nList):
        """ This looks at all entries in the provided coordinate lists
        and appends indices to the nList iff they are denoted as Well
        locations. """
        added = 0
        for x in xList:
            for y in yList:
                try:
                    curP = self.wellMatrix[x,y]
                    if curP == self.WELL: 
                        nList.append((x, y))
                        added += 1
                except:
                    continue
        return added
    
    def buildFeatureVector(self, rImg, gImg, well, k=8):
        """ This returns a k-tuple (8 by default; adjust with parameter k):
        [(C,I,D), ... , (C',I',D')]
        Where C is the colour of the nth neighbor, each I is the 
        average intensity of a neighbor to w (0 if considered "black"), and D
        is the Euclidean distance (in pixels) to that neighbor """
        neighbors = self.getKNearestNeighbors(well, k)
        vector = []
        for n in neighbors:
            nTup = []
            gIntensTup = self.getAvgIntensityForWell("green", n, gImg)
            rIntensTup = self.getAvgIntensityForWell("red", n, rImg)
            # If we're in conflict, ignore well
            if (gIntensTup[1] > 1 and rIntensTup[1] > 1) or (gIntensTup[1] <= 1 and rIntensTup[1] <= 1):
                nTup.extend(("black",0))
            elif gIntensTup[1] > 1:
                nTup.extend(("green", gIntensTup[0]))
            else:
                nTup.extend(("red", rIntensTup[0]))
            nTup.append(analyze.getEuclidDistance((well.x, well.y), n))
            vector.append(nTup)
        return vector
            
    def drawMapImage(self, outPath):
        RED = (255, 0, 0)
        GREEN = (0, 255, 0)
        print "Drawing well map image..."
        if len(self.redWells) == 0:
            print "Warning: no red wells. Make sure to build the talker list."
        if len(self.greenWells) == 0:
            print "Warning: no green wells. Make sure to build the green well list."
        else:
            print "Red well count: " + str(len(self.redWells))
            print "Green well count: " + str(len(self.greenWells))
        outImg = Image.new("RGB", self.sytosImg.size)
        for w in self.redWells:
            for xPrime in (0,1,2):
                for yPrime in (0,1,2):
                    try:
                        outImg.putpixel((w.x+xPrime, w.y+yPrime), RED)
                    except:
                        print "Pixel is out of bounds... ignoring"
                        continue
        for w in self.greenWells:
            for xPrime in (0,1,2):
                for yPrime in (0,1,2):
                    try:
                        outImg.putpixel((w.x+xPrime, w.y+yPrime), GREEN)
                    except:
                        print "Pixel is out of bounds... ignoring"
                        continue
        print "Finished. Writing image out to file..."
        outImg.save(outPath, "TIFF")
        print "Image map written out.\n"
        
    def getAvgIntensityForWell(self, colour, w, img):
        """ 
        Returns the average (pixel) intensity for the well whose
        upper left is specified by the w-tuple parameter. This assumes
        that img is grayscale! Average intensity is taken to be the
        average numerical amount above corresponding tNot value (i.e., difference from
        intensity at time=0). A scalar is returned that represents the average 
        intensity of the well.
        
        Requisite parameters:
        colour -- if this is "green" then the green threshold
                    scalar is used. it is assumed the img is measuring
                    intensity of green waves.
        w -- well coordinate for which average intensity will be returneed
                w represents upper left of 3x3 well.
        img -- pass in the green image if "green"; else pass in red img. 
                again, these are assumed to be grey scale intensity images.
        """
        total = 0
        numPixelsAboveThreshold = 0
        rVal = -1
        for xPrime in (0,1,2):
            for yPrime in (0,1,2):
                if colour == "green":
                   try:
                       p = img.getpixel((w.x+xPrime, w.y+yPrime)) 
                   except: 
                       # We're out of bounds here, just ignore
                       p = 0
                   total += p
                else:
                    # Then do it for reds.
                    try:
                        p = img.getpixel((w.x+xPrime, w.y+yPrime))  
                    except:
                        p = 0
                    if p>0:
                        total += p
        if colour == "green":
            return (total/9)-self.getTNotIntensityForWell(w)
        else:
            return (total/9)
          
                        

    